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Turbulent three-dimensional MHD dynamo model 
in spherical shells: Regular oscillations of the 

dipolar field 

By R. D. Simitev*, F. H. Busse§ and A. G. Kosovichev^, 



We report the results of three-dimensional numerical simulations of convection-driven 
dynamos in relatively thin rotating spherical shells that show a transition from an strong 
non-oscillatory dipolar magnetic field to a weaker regularly oscillating dipolar field. The 
transition is induced primarily by the effects a stress-free boundary condition. The vari- 
ation of the inner to outer radius ratio is found to have a less important effect. 



1. Introduction 

The Sun possesses a deep layer of intensely turbulent convection near its surface. 
The dynamical properties of this zone continue to elude basic physical explanation. The 
most important questions involve the solar differential rotation profile and the period- 
icity of the solar magnetic activity. The differential rotation is known from helioseis- 



mological studies (e.g. Schou et al., 19981: the angular velocity variation observed at 
the surface where the rotation is faster near the equator and slower near the poles, 
extends through the convection zone with little radial dependence. The basic feature 
of the solar activity is the existence of a periodic 22-year long cycle which is mainly 
manifested in the periodic appearance of sunspots accompanied by an oscillation of the 
antisymmetric (dipolar) parts of the solar magnetic field during which it reverses po- 



larity (e.g. Brandenburg & Subramanian, 20051. Both the differential rotation and the 



solar activity cycle are persistent global scale features of the Sun and it is difficult to 
understand how they emerge from the intensely turbulent fiow in the solar convection 
zone which includes scales ranging from granules (1 Mm in horizontal size), to super- 
granules (30 Mm), and possibly giant cells (over 200 Mm) and which shows similar 
contrasts in time scales. A variety of modeling approaches ranging from linear theory 
(e.g. Riidiger & HoUerbach, 2004) and weakly nonlinear theory ( |Busse, 1970[ ) to high- 



resolution three-dimensional global simulations of compressible convective dynamos in 
spherical shells ( |Brun fc Toomre, 2002[ |Browning, 2007[ ) have been employed but the 
precise mechanisms leading to the peculiar form of the profile of solar differential rota- 
tion, and the regularity of its magnetic activity remain unclear. 

In a recent study |Goudard &: Dormy, 2008 carried out three-dimensional numerical 



simulations of self-excited convective dynamos with the aim of exploring the dependence 
of solutions on the thickness of the spherical shell. They employed a model that has been 
used as a benchmark by the geodynamo community ( [Christensen et al., 2001"] ), which 
they modified by replacing the no-slip boundary conditions at the top of the spherical 
shell by stress-free ones. They found that a transition occurs from non-oscillating to 
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regularly oscillating dipolar dynamo solutions as the thickness of the spherical shell is 
decreased and concluded that this is a purely geometric effect. They discussed their 
finding in the light of the facts that the outer core, the layer where the non-oscillatory 
geodynamo is believed to be generated is a relatively thick layer in the case of the Earth, 
whereas the solar convective zone is a relatively thin. 

However, it is known that oscillating dipolar magnetic solutions can also be obtained 
in dynamos in thick spherical shells (Busse & Simitev 2006), and it remains unclear 
whether the transition found by Goudard & Dormy (2008) is the result of variation of 
the shell thickness or the introduction of a stress-free boundary condition. In this report 
we attempt to clarify this question and extend the results of Goudard & Dormy (2008) by 
exploring the effects of a different set of boundary conditions. We find that the transition 
to an oscillatory state is more likely due to the use of stress-free boundaries. 



2. Mathematical formulation of the problem and methods of solution 

We consider a spherical fluid shell rotating about a fixed vertical axis. We assume that 
a static state exists with the temperature distribution 

Ts^To + l:^Trir-^{l-ri)-'\ 

where r denotes the distance from the center of the spherical shell, -q denotes the ratio 
of inner to outer radius of the shell and d is its thickness, and AT is the temperature 
difference between the boundaries. The gravity field is given by 

g = —djr. 

In addition to d, the time d'^ /v, the temperature v'^ /jad'^, and the magnetic flux density 
v{^qY^'^ / d are used as scales for the dimensionless description of the problem where v 
denotes the kinematic viscosity of the fluid, k its thermal diffusivity, g its density and ^ 
is its magnetic permeability. The equations of motion for the velocity vector it, the heat 
equation for the deviation O from the static temperature distribution, and the equation 
of induction for the magnetic flux density B are thus given by 

dtU + u- Vu + Tk-Ku^ -Vtt + Gr + V^tt + B ■ VB, (2.1a) 

V-M = 0, (2.1&) 

V{dtQ + u ■ Ve) = (R?7r~3(l - ii)~'^)r ■ u + V^G, (2.1c) 

V-B = 0, (2. Id) 

V^B ^¥,^{dtB + U- VB ~ B - Vu), (2.1e) 

where dt denotes the partial derivative with respect to time t and where all terms in 
the equation of motion that can be written as gradients have been combined into Vtt. 
The Boussinesq approximation has been assumed in that the density g is regarded 
as constant except in the gravity term where its temperature dependence, given by 
a = —{dg/dT)/g =const, is taken into account. The Rayleigh numbers R, the Cori- 
olis number r, the Prandtl number P and the magnetic Prandtl number Pm are dcflned 

VK V K X 

where A is the magnetic diffusivity. Because the velocity field u as well as the magnetic 
flux density B are solenoidal vector fields, the general representation in terms of poloidal 
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and toroidal components can be used 

M = V X (Vi; X r) + Vw x r , (2.3a) 
S = V X (V/i X r) + Vg X r . (2.36) 

By multiplying the (curl)^ and the curl of equation (|2.1ap by r we obtain two equations 
for V and w 

[(V^ - dt)C2 + Td^]V^v + tQw - C2Q = -r • V X [V X (m • Vm - B • VB)], (2.4a) 
[(V^ - dt)C2 + Td^]w - tQv ^ r V X {u Vu - B ■ VJB), (2.46) 

where 9^ denotes the partial derivative with respect to the angle ip of a spherical system 
of coordinates r,9,ip and where the operators C2 and Q are defined by 

£2 = ~r^y^ + dr{r^dr), 
Q = rcoseV^ - {C2 + rdr)icos9dr - r'^ smede). 

The heat equation for the dimensionless deviation from the static temperature distri- 
bution can be written in the form 

V^e + {K'qr~^{\ - r])-'^)L2V = V{dt + u ■ V)e, (2.5) 

and the equations for h and g are obtained through the multiplication of equation (|2.1ep 
and of its curl by r 

V^/Iz/i = Pm[at/:2/i - r • V X (u X B)], (2.6a) 
¥^£23 = Pm[at/:25 - r • V X (V X (m X B))\. (2.66) 

In this report we discuss results obtained with two different sets of velocity boundary 
conditions, namely no-slip conditions given by 

V ~ drV ~ w = at 7' — ?'i = 77/(1 — 1]) and r = 7'o = 1/(1 — 77), (2.7) 

and mixed conditions, i.e. a combination of no-slip conditions at the inner spherical 
boundary and a stress-free condition at the outer boundary, given by 

V = drV = ui = at r = 7'i = 1]/ (1 — 77), (2.8) 
V = d^^v = dr{w/r) =0 at 7- = r^, = 1/(1 — 7/). 
As boundary conditions for the heat equation, we assume fixed temperatures 

e = at r r, = 77/(1 -77) and r = To = 1/(1 -77). (2.9) 

For the magnetic field, electrically insulating boundaries are assumed such that the 
poloidal function h must be matched to the function h^'^\ which describes the poten- 
tial fields outside the fluid shell 

g^h- h^"^ = drih - h^"'^) = at r = ri = 77/(1 - 77) and r = = 1/(1 - 77). (2.10) 

But computations for the case of an inner boundary with no-slip conditions and an 
electrical conductivity equal to that of the fluid have also been done. The numerical 
integration of equations (|2.4p . (|2.5p . and (|2.6p together with boundary conditions (|2.9p . 
(|2.10p and any of (|2.7p or (|2.8p . proceeds with the pseudo-spectral method as described 
by Tilgner (1999) which is based on an expansion of all dependent variables in spherical 
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harmonics for the 9, (/^-dependences, i.e. 

V = ^Vr'{,r,t)Pl^{cose)e^^{imtf] (2.11) 

and analogous expressions for the other variables, it;, 9, h and g. Here P™ denotes the 
associated Legendre functions. For the r-dependence expansions in Chebychev polyno- 
mials are used. For the computations to be reported in the following a minimum of 41 
collocation points in the radial direction and spherical harmonics up to the order 128 
have been used. 

3. Dependence of convection and dynamo solutions on the thickness of the 
spherical shell 

In order to clarify whether the transition from non-oscillatory to oscillatory dynamos 
reported by Goudard & Dormy (2008) is the exclusive result of variation in the radius 
ratio ry, we started with a case similar to the one described in their work and identical 
to the geodynamo benchmark simulation described in (Christensen et al. 2001). More 
precisely, we employ no-slip boundary conditions, as given in the preceding section, both 
at the inner and outer surface of the shell. We then simulate a sequence of dynamo 
solutions where we increase the value of 77 while keeping all other parameter values fixed. 

Global quantities that characterize the basic properties of dynamo solutions include 
the magnetic energy density that can be decomposed into several components, namely 
poloidal and toroidal energy densities, each of which can be further decomposed into 
mean and fluctuating parts. The corresponding definitions are given by the expressions 

Vx (VT^xr) = Vgxr 

Vx (V/ixr) p), A/, = i(| Vgxr p), 

where (•) indicates the average over the fluid shell and h refers to the axisymmetric 
component of h. h is deflned hy h = h -- h. Similarly, kinetic energy densities Ep, Et, 
Ep, and Et can be defined with analogous expressions where v and w replace h and g. In 
addition, the energy densities can be divided into those of fields that arc antisymmetric 
(axial dipole symmetry) and those that are symmetric (axial quadrupole symmetry) with 
respect to the equatorial plane. The dipole (quadrupole) fields are described by spherical 
harmonic Y™ with odd (even) I + m. Other global quantities of interest are the helicity 
defined as 

He = ((V X u) ■ u), 
and the cross-helicity which is given by 

XHe = {u ■ B). 

In Fig. [1] these first-order characteristics of dynamo solutions are plotted as a function 
of the radius ratio rj of the spherical shell, whereas in Fig. [5] snapshots of typical spa- 
tial solution structures are presented. We remark that a value of ry = 0.35 is typically 
assumed to be relevant in the case of the Earth, where the inert inner core extends to 
less than 40% of the core radius; a value of rj ~ 0.7 is appropriate for the Sun, where the 
radiative zone fills about 70% of the solar radius. The well-known benchmark solution 
of (Christensen et al., 2001) has 77 = 0.35 and represents a strong dipole. In Fig. [1] this 
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Figure 1. The averaged magnetic, M, and kinetic, E, energy density components, and the 
helicity. He, as a function of the radius ratio rj for fixed values of the other parameters: P = 1, 
R = 10^, r = 2000, Pm = 5 and no-slip velocity boundary conditions. The components Xp, Xt, 
Xp, and Xt are shown by crosses, squares, triangles and circles, respectively, where X = M,E. 
The energy density Ep, has not been included because it is nearly two orders of magnitude 
smaller than Et- The values of helicity are indicated by stars. 



is evident from the fact that the mean poloidal dipolar magnetic energy Mp makes the 
dominant contribution. This solution is known to be quasi-stationary in that the energy 
densities remain constant in time, and the only time dependence is the drift of the con- 
vection pattern in the azimuthal direction. The solutions at larger values of the radius 
ratio 7] exhibit a similar behavior, namely they remain strong non-oscillatory dipoles. 
All components of the kinetic energy, as well as the helicity, grow by roughly an order 
of magnitude as -q increases from 0.35 to 0.65. The magnetic energy components grow 
only slightly despite the increasingly vigorous convective flow, and at 77 = 0.65 dynamo 
action is lost indicating that as the shell thickness decreases it is more difficult to sus- 
tain magnetic field generation. The relative contributions of the various magnetic energy 
components remain unchanged. An inspection of Fig. [2] indicates that the main effects 
of increasing rj are the growth of the wavenumbcr of convection and the migration to 
lower latitudes of both, convective motions and magnetic fields, while the polar regions 
gradually become less active. To ensure that the described behavior is not transient we 
have continued the simulations for more than 20 magnetic diffusion times in each case. 



4. An attempt to model the solar magnetic cycle 

It appears from the results presented in the preceding section that variations of the 
shell thickness alone are not sufficient to induce oscillations in dipolar dynamo solutions. 
It is known that dipolar oscillations may be excited in a variety of other situations 
(Busse & Simitev 2006). An oscillating dipolar solution can be found when a region 
in the parameter space is approached where the quadrupolar or the higher-multipole 
components of the magnetic field arc not negligible. As these components are typically 
oscillatory, an oscillation in the dipolar parts is also excited. A region of multipolar 
dynamos may be approached, for instance, by reducing the value of the magnetic Prandtl 
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Figure 2. Snapshots of selected components of the dynamo solutions with P = 1, R = 10^, 
r = 2000, Pm = 5, and no-slip velocity boundary conditions for varying radius ratio 77. Each 
circle of the first column shows lines of constant in the left half and of r sin Odgh = const, in 
the right half. The second column shows meridional lines of constant in the left half and of 
r sin Odgv in the right half. The third column shows meridional lines of constant values of the 
azimuthally averaged cross-helicity in the left half and helicity in the right half. Each plot in 
the last column shows poloidal streamlines rd,pV = const, in the equatorial plane. The value of 
r) = 0.35, 0.5, 0.6 in the first, second, and third row, respectively. 

number Pm, or by increasing the value of the Rayleigh number R or of the rotation 
parameter t. In all of these situations one finds that the increase of the differential 
rotation plays a crucial role in exciting the oscillations. The differential rotation may, 
of course, be enhanced much more readily by imposing a stress-free boundary condition 
for the velocity field. Thus, the suggestion of Goudard & Dormy (2008) to replace the 
no-slip condition on the outer spherical boundary by a stress-free one is vifell-justified. As 
these authors remark, the choice of mixed velocity boundary conditions may be argued 
also on physical grounds: the no-slip condition at the inner boundary mimics the solar 
tachoclinc (Tobias 2005) and the stress-free condition at the outer boundary mimics the 
solar photosphere. 

In Fig. [3] a dynamo simulation with mixed velocity boundary conditions is presented. 
After an interval of about 5 viscous diffusion times an abrupt transition in the nature of 
the dynamo solution occurs. The most notable signature of the transition is the decrease 
of the magnetic energy density by more than an order of magnitude. At the same time 
the kinetic energy increases tvifofold, mostly due to the increase in differential rotation. 
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Figure 3. Time series of a dynamo solution in the case P = 1, R = 10^ r = 2000, Pm = 4.5, 
?7 = 0.65 and mixed velocity boundary conditions. The first row shows dipolar magnetic energy 
densities. The second row shows kinetic energy densities. The last row shows values of the mean 
dipolar, H^, and the mean quadrupolar, H2, components in the spherical harmonic expansion 
of the magnetic field. The components Xp Xt, Xp, and Xt are indicated by circles, squares, 
plus signs and crosses attached to the respective curves, and X stands for either M or E. The 
coefficient Hi is represented by a solid line and H2 by a dashed line. The time series in the right 
column represent enlarged sections of the ones in the left column. 



This is not surprising as the increase is promoted by two effects: first, by the removal of 
the no-slip condition and second, by the decrease of magnetic field vifhich, now becomes 
too weak to effectively act as a brake on differential rotation. This latter effect has been 
reported in numerous other cases (Bussc 2002, Simitcv & Bussc 2005). 

More subtly, the transition is observed as a change in the relative contributions of 
magnetic energy components. The mean poloidal dipolar magnetic energy Mp is no longer 
the dominant component and is overcome by the fluctuating parts of the magnetic energy. 
In this respect the transition is similar to the transition between Mean Dipolar (MD) 
and Fluctuating Dipolar (FD) dynamos recently reported by Simitev & Busse (2009). 

The temporal behavior of the dynamo changes dramatically after the transition as 
well. In particular, it exhibits nearly periodic dipolar oscillations. These are evident in 
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Figure 4. Half a period of dipolar oscillations of the dynamo shown in Fig. [S] Each circle of the 
first column shows lines of constant in the left half and of rsmOdeh = const, in the right 
half. The second column shows meridional lines of constant Br at r = ro. The third column 
shows meridional lines of constant u,p in the left half and of rsin69eTJ in the right half. Plots 
follow from top to bottom with a time step At — 0.021665. 



the right-hand side of Fig. [3] and even better in the oscillations of the spherical harmonic 
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Figure 5. A butterfly diagram for the case shown in Figs. [3]and|4l i.e. the value of rdeg at 
r — To — 0.05 as a function of time and latitude. 

coefficient describing the contribution of the axial dipole in the solution. It is also 
remarkable that the quadrupolar coefficient H2 remains much weaker than in correspond- 
ing oscillations in thick shell dynamos. A sequence of snapshots equidistant in time is 
shown in Fig. U in order to visualize half a period of oscillation. Except for the sign the 
plots shown in the last row closely resemble the ones in the first row confirming thus the 
near perfect periodicity of the oscillation. A polarity reversal of the magnetic field occurs 
during the cycle and is evident both in Fig. |4] and in the change of sign of in Fig. [3l 
This cyclic behavior is sustained for over 25 viscous diffusion times although it appears 
from the right-hand side of Fig. [3] that a weak modulation with a second frequency is 
also present. Finally, we remark that the oscillations appear to be well-described by the 
Parker wave linear analysis presented in Busse & Simitev (2006), in that the period of 
oscillation T = 0.12996 is well approximated by the expression 



which in this case yields 0.11036. 

In order to stress the similarity with the solar magnetic activity we present in Fig. [5] 
a butterfly diagram of the simulation, where the cyclic behavior and the propagation of 
magnetic features from high latitudes toward the equator is also visible. 

5. Discussion and outlook 

We have confirmed the finding of Goudard & Dormy (2008) of a sharp transition 
from a dynamo dominated by a non-oscillatory dipolar magnetic field to a nearly perfect 
oscillatory dynamo with a weaker dipole. We have related this transition to corresponding 
transitions from dynamos with a strong nearly steady dipole to oscillatory dynamos in 
thick shells. The introduction of a stress-free boundary condition at the outer spherical 
surface rather than the reduction of the shell thickness appears to be essential for the 
transition. It is of interest to note that the oscillations of thin shell dynamos are much 
more purely dipolar than those found in thick shells. 

It is tempting to compare the oscillatory dipolar dynamo of section 3 to the magnetic 
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activity of the Sun. A notable difference with the Sun is the profile of differential rotation 
of the solution which has the form of a nearly constant angular velocity on cylinders. 
Although this profile, when mapped to the surface r = Tq, reproduces qualitatively the 
rotation of the solar surface, the differences in radial direction arc significant. 

This is not surprising in view of the strong change of density throughout the solar 
convection zone caused by compressibility which has not been taken into account in 
our simple model. In addition the scales of turbulence are far from being resolved. This 
emphasizes the need to explore the parameter space of the problem in further detail. 

We gratefully acknowledge support from CTR and NASA which made our visit to 
Stanford possible. The numerical calculations reported in this paper were carried out us- 
ing the computer resources of the School of Mathematics and Statistics of the University 
of Glasgow, and the UK MHD supercomputer at the University of St. Andrews, UK. 
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